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ABSTRACT 

Small perturbations in spherical and thin disk stellar clusters surrounding massive a black hole are 
studied. Due to the black hole, stars with sufficiently low angular momentum escape from the system 
through the loss cone. We show that stability properties of spherical clusters crucially depend on 
whether the distribution of stars is monotonic or non-monotonic in angular momentum. It turns 
out that only non-monotonic distributions can be unstable. At the same time the instability in disk 
clusters is possible for both types of distributions. 
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1 INTRODUCTION 

The study of the gravitational loss-cone instability, a far 
analog of the plasma cone instability, has begun with the 
work of V. Polyachenko (1991), in which a simplest analyt- 
ical model of thin disk stellar cluster has been treated. The 
interest to the problem of stability of stellar clusters has 
been revived recently by detailed investigation by Tremaine 
(2005) and Polyachenko, Polyachenko, Shukhman, (2007; 
henceforth, Paper I) of low mass clusters around massive 
black holes. The both papers have considered stability of 
small amplitude perturbations of stellar clusters of disk-like 
and spherical geometry. 

Tremaine (2005) has shown using Goodman's (1988) 
criterion that thin disks with symmetric DFs over angu- 
lar momentum and empty loss cone are generally unstable. 
By contrast, analyzing perturbations with spherical num- 
bers I — 1 and I = 2, he deduced that spherical clusters with 
monotonically increasing DF of angular momentum should 
be generally stable. 

Later we demonstrated (see Paper I) that spherical sys- 
tems with non-monotonic distributions may be unstable for 
sufficiently small-scale perturbations / ^ 3, while the har- 
monics I = 1, 2 are always stable. For the sake of conve- 
nience, we have used two assumptions. The first one is that 
the Keplerian potential of the massive black hole dominates 
over a self-gravitating potential of the stellar cluster (which 
does not mean that one can neglect the latter). Then the 
characteristic time of system evolution is of the order of the 
orbit precessing time, which is slow, compared to typical 
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dynamical (free fall) time. Since a star makes many revo- 
lutions in its almost unaltered orbit, we can regard it as to 
be "smeared out" along the orbit in accordance with passing 
time, and study evolution of systems made of these extended 
objects. 

The second assumption is a so called spoke approxima- 
tion, in which a system consists of near-radial orbits only. 
This approximation was earlier suggested by one of the au- 
thors (Polyachenko 1989, 1991). The spoke approximation 
reduces the problem to a study of rather simple analytical 
characteristic equations controlling small perturbations of 
stellar clusters. 

There are two questions that naturally arise in this con- 
text. First: Does the instability remain when abandoning the 
assumption of strong radial elongation of orbits? Second: 
Does the instability occur in spheres with monotonically in- 
creasing distributions in angular momentum if one consider 
smaller-scale perturbations with I ^ 3? The aim of the paper 
is to provide answers to these questions. 

To achieve the task we use semi-analytical approach 
based on analysis of integral equations for slow modes elab- 
orated recently in Polyachenko (2004, 2005) for thin disks, 
and in Paper I for spherical geometry. Following Paper I, we 
shall restrict ourselves to studying monoenergetic models 
with DFs in the form 

F{E,L) = A5{E-E )f{L). (1.1) 

The models specified by function f(L) are suitable for study- 
ing the effects of angular momentum distribution on gravi- 
tational loss-cone instability. On the other hand, the Dirac 
5-function permits one to reduce the integral equations for 
slow modes to one-dimensional integral equations, and to 
advance substantially in analytical calculations. 
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Several arguments can be brought in favour of our sim- 
plified approach. First of all, the Lynden-Bell derivative (see 
Paper I, eq. 4.7) of the DF with respect to angular momen- 
tum L, keeping J = L + I\ constant (here 7i is the radial 
action) in the limit where the slow mode approximation is 
applicable, can be replaced by a derivative, keeping energy 
E constant: 



dF 
~dL 



dF 
dE 



+ 



dF 
~dL 



dF 



because Q pr is small. Thus, the derivative over energy is 
not included into the slow integral equation, and one can 
loosely say, that dependence on energy is only parametric. 
Another argument is that the results of independent study 
by Tremaine (2005), who used a non-monoenergetic DF, are 
in agreement with our conclusions. 

Section 2 is devoted to spheres, Section 3 - to thin disks 
with symmetric DFs. The sections are organized alike. In the 
beginning we derive integral equations for initial distribution 
functions in the form JTTT}. Then follow analytical and nu- 
merical investigations of these equations. We demonstrate 
that by contrast to the case of near-Keplerian sphere, the 
loss-cone instability in disks takes place even for the mono- 
tonic DF, df /d\L\ > 0, provided the precession is retrograde 
and the loss cone is empty: /(0) = 0. Sec. 2 is complimented 
by stability analysis of models with circular orbits, which of 
course doesn't belong to the class of monoenergetic models 
of <HHI) type. 

In the last, Section 4, we discuss the results and some 
perspectives of further studies. 



2 SPHERICAL SYSTEMS 

2.1 Integral equation for slow modes in 
monoenergetic models 

For the near-Keplerian systems, the slow integral equation, 
which has been derived in our Paper I (see there Eq. (4.8)), 
is neatly suited. In contrast to Paper I, we shall not assume 
here strong elongation of orbits, i.e. we shall go beyond the 
spoke approximation. 

Since energies of all stars are identical, the unper- 
turbed DF depends on one variable only. It is convenient to 
use a dimensionless angular momentum a = L/L c i IC (Eo), 
where L c i IC is the angular momentum on circular orbits: 
L c irc{E ) = GM C /(2\E \) 1/2 , M c is the central point mass, 
G is the gravitational constant. The frequency of stellar ra- 
dial oscillations (li(Eo) = (2\E \) 3/2 /(GM C ), and the radius 
of the system R(Eo) = GM c /\Eq\ are independent of the 
angular momentum. For shorthand notations, we shall omit 
the argument Eq. 

The normalization constant A is taken so that a mass of 
the spherical system surrounding the central mass is equal 
to Mg (we assume the ratio e = Mg/M c to be small: e<l): 



M G = / FdT = 2 (2tt) 



dE 



L ciTC 



LdL F(E,L). 



If one defines the normalization of the dimensionless DF over 
angular momentum / (see (1.1)) as L daaf(a) = 1, then 



normalization factor A in (1.1) is 



A = 



167T 3 L2, 



(2.1) 



It allows to represent the kernel of the integral equation 
(formula (4.8) in the Paper I) in the form 

P U,(E,L;E',L') = ^!i|±!) C,/C« ,(a,a'), 

where I is the index of the spherical harmonic, Ci = 
J °° dz z _1 [«/(;_)_ x -)/2 (z) Ji/2{z)} 2 and J v (z) is the Bessel 
function!]] The functions , satisfy to the condition 

JC^ , (0, 0) = 1; their explicit form is given later. Then sub- 
stitution of the DF in the form (II. ip leads to the following 
integral equation: 

i 

4>s{a) = 2ftieC; ]T s' 2 Df' x 

s '- s min 

1 

X ] u;*-s>*nUa>) K ' " (Q '°° Aa)da ' (2 - 2) 
o 

where <j) s (a) is the Fourier harmonics of the radial part of the 
perturbed potential (for more detail, see Paper I), fi pr (a) is 
the orbital precession rate, s m i n = 1 for odd I, and s m i n = 2 
for even I. The coefficients D are calculated by the formula 

_1 (i + s)!(I-8)! 

2 2l 



Df 



\l — s\ even, 



\l-s\ odd. 

(2.3) 

Recall that Eq. (|2.2|) is written in a noninertial reference 
frame centered on the mass M c . Then, additional indirect 
potential arising from the acceleration of the frame should 
be considered (see, e.g., Tremain 2005) 

g*(r,t) = GtJ r' 5p ^' 3 '^ dV, Sp = Jsfdv, 

where Sf is the perturbation to the background DF. Tremain 
(2005) argued that for the secular perturbations, this indi- 
rect potential must be omitted. Indeed, in studying secu- 
lar evolution one should consider perturbations 5f averaged 
over Keplerian orbits. In this case the perturbed density is a 
superposition of contributions of individual orbits, averaged 
over their periods. A special feature of a Keplerian orbit is 
that the average force from this orbit acting to the mate- 
rial point located in a focus of the ellipse is equal to zero. 
One must be careful however, since the perturbation is not 
well defined for orbits with low angular momenta. Below we 
shall consider systems with "small amount" of stars with low 
angular momenta only (see also discussion in Sec. 2.2.1). 

By changing the unknown function 

[w 2 — s 2 fipr(a)] ifi s (a) = </> s (a) 



1 For 1 = 1 the coefficient C\ can be calculated analytically: 
Ci = 4/3-7T 2 rj 0.135. Numerical calculations show decreasing C; 
with increasing the mode number I: Ci = 0.063, C3 = 0.0373, 
d = 0.025, C* 5 = 0.018 , and so on. 
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Eq. \Y2:2\\ can be reduced to the linear eigenvalue problem 

l 

[lu 2 - S 2 nl I {a)]i f , a (a) = 2n 1 £Ci ^ s' 2 D? x 

s ' — s min 

1 

x |n p( ( Q >'^M/C ( s ' ) s ,(a, a > J /(aVa' (2-4) 
o 

For almost radial orbits, when a < 1 or eccentricity 
e = \/l — a 2 ~ 1, the precession rate is 



fipr(a) = 



[l + 0(a 2 )]. 



(2.5) 



For orbits with smaller eccentricity, the numerical coefficient 
preceding the small parameter e Q\ is somewhat greater than 
2/tt 2 . Since one suggests that the characteristic frequencies 
of the problem under consideration are of the order of typical 
precession velocities, lu ~ fi pr ~ efli, it is convenient to 
change to the dimensionless frequencies, measured in the 
natural "slow" frequency: 

w= ^n? " (Q) = -^r- (2 - 6) 

For the spherical systems, the precession is retrograde (see 
Tremaine 2005, or Paper I), so v(a) > 0. Then the dimen- 
sionless integral equation becomes 



2 (a)]^(«) = -2C ; J2 S ' 2D 



i x 



v{a!) a' K. (, s ] s , (a, «')*>.' («') da', (2.7) 



To obtain the eigenfrequency spectrum for a model it is 
necessary to compute preliminarily the kernels K ,(a,a') 
(universal for all models), and the precession rate profile 
v(a) for the given model. The integration over Keplerian 
orbits is most conveniently expressed using the variable r, 
which is connected with the current radius r and the true 
anomaly C, of a staiQ as follows: 

r = i-R(l-ecosr), cosC= cosr ~ e . ( 2 .8) 
^ 1 — e cos r 

Then after some transformations, the kernel A^ ! ' , can be 

' s, s 

reduced to the form 



(21 + l)7r 2 Ci 



dr r cos(s^) X 



77 

J dr' r' cos{s'(')^i{r,r'), (2.9) 



is used. 

The expression for the precession rate can be obtained 
by transformation of expression (4.2) of Tremaine (2005) 
(see also Paper I): 

„(«) = -£- M0( cos i- e )<fr (2 . 10) 

47re J r 2 
o 

and !/(l) = — i/rp(i-R), where the density p(r) is defined 
by (232}. 

For the monoenergetic models, the minimal and maxi- 
mal radii are i? min = | -R(l-eDiax), -Rmax = | -R(l+e max ), 
where e max = (1 — h 2 ) 1 ^ 2 , and ft is the minimal dimension- 
less angular momentum corresponding to the boundary of 
the loss cone. 

The function p(r) is a ratio of the mass of a spherical 
system inside the sphere of radius r to the total mass Mg: 

T 

M(r) = n§^' ^°«= 47r / r' 2 p(r')dr', (2.11) 

fl mia 

Mo = -Mo^max), and the density is calculated by the for- 
mula 



p(r) = 



2a da /(a) 



(2.12) 



where a max = y^4 (r/i?)(l — r/R) . From here on we shall 
assume R = 1. 

Using (25} and (|2. 10[) — (j2. 12[) one can transform the 
expression for the scaled precession rate v(a) to a more uni- 
versal form: 

i 

= 5-1-3 fda'a'f(a')Q(a,a'), (2.13) 
o 

where the kernel Q(a,a') doesn't depend on a DF and 
equals to 



f * W-Z$Z-?r (2 - I4) 

Pmin 

with p m in = max (r m i n , r' min ), p max = min (r max , r^ax)- 
Here r min = \ (1 - e), r max = § (1 + e), r^ in = § (1 - 
e'), C ax = i(l + e '),ande=(l-a 2 ) 1 / 2 , e' = (1-a' 2 ) 1 / 2 . 
For the near radial orbits Q(0,0) = 4, so that one obtains 
the above mentioned result (|2.5[) : f w (2/7r 2 )a. 



where r' and £' specify the position of a star on the orbit 
with the eccentricity e corresponding to the variable t , and 
the notation 

■ / /\i 

Ti(r, r) 



max(r, r')' +1 

2 True anomaly is the angle between directions to the star and 
to the pericenter. 



2.2 Analytical results 

2.2.1 Exact solution with zero frequency for the lopsided 
mode (1 = 1) 

Tremaine (2005) has noted that for an arbitrary distribu- 
tion F(E, L) with empty loss cone, F(E, L = 0) = 0, 
a zero frequency lopsided mode I = 1 must exist. The 
mode corresponds to a non-trivial perturbation arising un- 
der shift of the spherical system as a whole relative to the 
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central mass. The perturbed potential in such a mode is 

<5$(r, 60 = — £cos# — ; — , where £ is the displacement. In 
dr 

terms of the function <f) s =i(a), this perturbation has a form 



or in terms of the function ipi (a) from (|2.7[) 
(pi(a) 



(2.15) 



(2.16) 



One can check that (|2.15[) and ui = provided the con- 
dition 

f(a = 0) = 



is a solution of (|2.2|l or (|2,7|l for 1 = 1, taking into account 
the expressions (I2.13P and (|2.14|l . written in the form 



Q{a,a ) = -i^ -J^- I dr y/ (r - r m i n )(r max - r) x 
a oa 



is p(0) = 0. In what follows we suppose this condition to 
be fulfilled. The conditon is not equivalent to the condition 
f(a = 0) = 0, imposed to the DF for the existence of such a 
solution of our governing integral equation, but it is equiv- 
alent to the stronger condition: f(a = 0) = /'(a = 0) = 0. 
Indeed, it is easy to show that if f(a) oc a 3 for small a then 
p(r) oc 7-( s_1 )/ 2 for small r. So, the condition s > 1 must be 
fulfilled. 

By contrast, in the disk systems the analogous m = 1 
zero mode does not exist, because there is no analog of the 
Newton's first theorem. 

The very existence of zero modes is crucial for stabil- 
ity analysis of spherical clusters with monotonic distribu- 
tions. Indeed, the role of destabilizing contribution of the 
second term in the right side of (|2.7p falls off with increasing 
the number I. So, it is expected that the most "dangerous" 
modes correspond to the lowest values of I. But it turns out 
that I — 1 mode is neutrally stable, and the next dangerous 
mode I = 2 is stable. Note, however, that such a reasoning 
is not valid for systems with non-monotonic distributions. 



and also the expression for the kernel K,^{ci,a' 



K$(a,a') = — - / dr\/(r- r min )(r max - r) x 



The lopsided solution with zero frequency is specific for 
spherical systems. At the first glance, it defies common sense 
to argue that the stationary mode in which the center of 
mass of a spherical system does not coincide with the BH 
is physical. Indeed, it seems that movement (oscillations) of 
the stellar cluster and the BH relative to the common center 
of mass must occur. However, it does not occur. 

It is, by all means, clear in the case of the empty loss 
cone of finite size, h > (here ft is a minimal value of di- 
mensionless angular momentum a, for which f(a) > 0). 
Indeed, let us consider the spherically symmetric cluster. 
Since the loss cone is finite, there is a spherical empty cav- 
ity of finite radius in the centre of sphere. According to the 
Newton's first theorem (Binney & Tremaine 1987), in this 
cavity the BH does not experience a net gravitational force 
from the cluster. Thus, if the BH is initially deposited at 
some arbitrary point within the cavity, it would remain at 
this position (and hence, acceleration of stellar cluster due 
to non-coincidence of centers of mass does not occur) . 

In the case when h — the situation is not so obvious, 
but the net force acting to the BH from the shifted spherical 
system can be zero as well. In order to assure this, one should 
write down the indirect potential taking into account the 
expression for perturbed density in zero lopsided mode Sp — 
— £ cos Qdpjdr: 



$*(r,t) = ~2n£Gr cosfl / dr 



4tt 



j dp(r') 
dr' 



Gr cos0£p(O), 



'■9' sm6' d0' 



2.2.2 The stable mode in systems with near-radial orbits 

By analyzing (|2.7|l . it is easy to find one more analytical 
solution with the frequency Gj = 0(1) at arbitrary values of 
I, for the models with highly elongated orbits. First of all 
we note that the frequency of this stable mode corresponds 
to the resonance on the tail of a narrow distribution, and so 
it decays exponentially slowly. In this way the mode differs 
from the unstable modes of interest which have a resonance 
in a region where the distribution is localized, i.e. at a < olt\ 
so they have characteristic frequencies and growth rates of 
the order of 0(o.t). 

After setting u> ~ 1 ^> «t in (12. 7p . omitting the second 
summand in l.h.s., turning to the spoke approximation, and 



taking into account the equality 
one finds 



E s 2 D? = U(l + l), 



(2.17) 



Hence, the condition for omitting of the indirect potential 



It is essential that this high-frequency mode is independent 
of details of the DF. Note also that in the systems with pro- 
grade precession, this mode describes the well-known radial 
orbit instability (instead of the neutral oscillations). 



2.2.3 The Variational principle 

Using l|2.7fl . one can prove two important statements: 

(i) For spherical system models with monotonic distribu- 
tions /(a), the eigenfrequency squared, uj 2 , must be a real 
number. This means that either the instability is absent at 
all, or aperiodic instability with RetD = occurs. 

(ii) Rotating (or oscillating) unstable modes may appear 
only in models with non-monotonic distributions. 
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Let us write Eq. (|2,7[1 in the form 



- 2 



(fis(a) = s 2 v 2 (a)ip 3 (a) — 2Ci^2 s ' 2 D i x 

s' — S m j n 

1 

x j g{a')K.^ sl {a,a')vs>{a')d(J, (2.18) 



where g(a) = v(a)a df(a)/da. We multiply both parts of 
Eq. (|2. 18fl by s 2 Dfipl (a), sum the result over s (asterisk 
means the complex conjugation), and integrate over a with 
the weight g(a). Then we obtain 



- 2 



/ g(a)da 22 s 2 Di \tp s (a)\ 2 = 

o s=s mi „ 

r 

= f 2 (a)g(a)da ]P s4d i \fs(a)\ 2 - 
- 2d J da J da J2 ( ss ') 2 A" Df : 

s — s min s'— = ■ 



x /C«,(a,a') fo(a)rf(a)] [ 5 (a') ^(a')]- (2.19) 

The reality of the coefficients of a) 2 in the l.h.s. of (|2.19[) and 
the first term in the r.h.s. is evident. With the help of (|2.9p . 
one can show that the kernel in Eq. (12 . 1 9 p has the following 
property of symmetry: 

4 ! > s ,(aV) (2.20) 

So, one can readily see that the second term in the r.h.s 
is real also. Consequently, taking the imaginary part of Eq. 
(|2.19|) . one obtains 



2 f ' 2 

Im(oi 2 ) / g(a)da ^ s 2 Dt 



\Vs(a)\ 



(2.21) 



^From the last equality, the statements formulated 
above follow immediately. If the function g(a) (or, equiv- 
alently, df(a)/da) is constant-sign, then the integral should 
be non-zero, and so Im(u> 2 ) = 0. In contrast, when Im(tj 2 ) 7^ 
0, the integral must be equal to zero. Consequently, the 
function g(a) should change its sign, i.e. DF f(a) is non- 
monotonic. 

Let us explain the term variational principle used in 
the title of this subsection. Consider a dynamic equation in 
the form d 2 £/dt 2 = — lu 2 £ — —K£. Provided that "elasticity 
operator" K is Hermitian, the dynamic equation may be 
obtained from the conditions 8{lu 2 )/8^ = and 6(lu 2 )/S£* — 
0. Here <5£ and <5f * are considered formally as independent 
variations while the functional uj 2 is 

j _ !Cm)w{a)da 
J \£\ 2 w(a) da 

(w(a) is a nonnegative weight function). In such a case it 
is used to speak about the variational (or energy) princi- 
ple (see, e.g., review by Kadomtsev 1966 on MHD-stability 
of plasma). From the other hand it is easy to see that if 
/ \£,\ 2 w{a) da^O for any nontrivial £ then uj 2 is real. Thus 
usually (as is the case in MHD-stability theory of plasma 
where K is Hermitian and w > 0) the notions "variational 



principle" and reality of uj 2 are identical. However, in our 
case the condition f \£\ 2 w(a) da 7^ is not satisfied for any 
nontrivial £ unless the DF is monotonic. Assuming that this 
condition is fulfilled, following the tradition that originates 
from plasma physics we speak that the variational principle 
takes place. Only in this case the dynamical equation can 
be interpreted mechanically, in terms of elastic forces. 

Evidently, the condition (|2.21|l is a serious obstacle to 
obtain unstable rotating modes. So, one might want to get 
round this obstacle. For instance, if we slightly change the 
initial monotonically increasing DF in a narrow region near 
a — 1, to vanish quickly but smoothly, then a modified sys- 
tem would be practically indistinguishable from the initial 
one. But then the variational principle breaks down. The 
question can be formulated as follows: May the discontinu- 
ous vanishing of f(a) at a = 1 be considered as the violation 
of monotony for the DF? 

Importance of this question is known since stability 
study of stellar systems with isotropic DFs, F = F(E) 
(Antonov, 1960, 1962). The variational principle there re- 
quired a DF to be decreasing function of energy E, F'(E) < 
0, everywhere. The systems with F'(E) > need separate 
examinations, that was carried out in some cases (see, e.g., 
Antonov, 1971, Kalnajs, 1972, Polyachenko and Shukhman, 
1972, 1973, Fridman and Polyachenko, 1984). An essential 
difference of the second type of DFs is in jumps to zero at 
the phase space boundary E = -Ebound- In fact, we have 
in this case an interval degenerated into the single point 
E — -Ebound where F(E) is decreasing. 

We checked numerically a possibility of the instability 
development connected with the maximum on the edge of 
the distribution function's domain. For this purpose, num- 
ber of models smoothed near a — 1 were computed. The 
computations showed no sign of instability, in contrast to 
isotropic distributions, F — F(E). The reason for the tol- 
erance of our present models is in fact that the kernels K, 
of integral operators in (|2.18|l vanish for the circular orbit 
a — 1, so details of the initial distribution near circular or- 
bits cannot affect much solutions of the integral equation 
(pT8)l . 

Roles of different terms in Eq. (|2. 18f) can be easily un- 
derstood. When dF/dL > 0, the first term of the right side 
in Eq. (I2.18P favors stabilization, while the second term gives 
destabilization (taking into account that the operator in- 
volved into this construction is self-adjoint and positively 
denned). In principle, this destabilizing effect could lead to 
instability. However, this is not the case because the sta- 
bilizing contribution exceeds destabilizing one in all cases 
considered both by Tremaine (2005) , and in the present pa- 
per (see the following sections). 



2.3 Unstable models 

Instability boundaries in terms of the angular momentum 
dispersion qt < (ar) c found in Paper I for the monoener- 
getic DF with 



er \ N ( a 2 1 . ■. 

JW ^TT TT ex P(-" / q t), 



(2.22) 



(N is the normalization constant, ar is the dimension- 
less angular momentum dispersion, n is the real number) 
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have a qualitative character only: formally, these bound- 
aries lie outside the validity of the spoke approximation, 
since (qt)c ~ 1. Obtaining such critical dispersions means 
only that the spoke models, in which at <C 1 by definition, 
are certainly unstable. So the quantitative determination oi 
these boundaries with help of the exact integral equation 
is required. The power -exp model (|2.22p is studied in Sec. 
2.3.1. 

In Sec. 2.3.2 we study a simplest Heaviside model con- 
sisting of two steps (at a = h\ and at a — h 2 ) (both in 
the spoke approximation framework and using exact inte- 
gral equation): 

/(«) 



h\ - ht 



In(a7^)exp(-a7a4), 



z n exp( — z) dz, z=—~-. 



o 



p(r)=N 



I 



' dz 



2 [H(a-hi)-H(a-h 2 )\, hi<h 2 <l 

(2.23) 



(H(a) denotes the Heaviside function). Finally, in Sec. 2.3.3 
we consider the log - exp model with DF 

N 



(2.24) 



for a ^ h, and f(a) = for a < h, i.e., with the empty loss 
cone (N is the normalization constant). 



2.3.1 The power - exp model 

Following Paper I, here we consider the stability of models 
(|2.22[) with n — 2 and n = 3 relative to the spherical har- 
monic I — 3. We remind that at the limit ot >C 1, both these 
models were unstable (the stability boundaries obtained us- 
ing spoke approximation were («t)c = 0.193 for n — 2, and 
(a T ) c = 0.283 for n = 3). 

For distribution (|2.22|l one finds 



Particularly, in the case «t <C 1, the normalization constant 
is N = 2/(n !) . From (|2J2|| , we obtain 



Further calculations of the density (|2.12p and precession rate 
(|2 . 1 3|) profiles should be evaluated numerically. 

Solutions of the integral equation (|2.7|l for n — 2 and 
n = 3 show similar behavior. At small values of aj, each 
model has one unstable mode. With increasing the dimen- 
sionless angular momentum dispersion cct, the growth rate 
of the instability decreases, and then vanishes at some criti- 
cal value (qt) c : for the model n — 2, (qt)! 2 ' — 0.301, for the 
model n — 3, (qt)! 3 ' — 0.311 (see Fig.[T|. We conclude that 
the spoke approximation in this case is qualitatively correct, 
but quantitatively poor. The instability becomes saturated 
at some critical value (qt)c, while the discrepancy between 
exact and approximate values of {olt)c for both models are 
not small. 

Apart from the unstable mode, the spectrum of each 
model has a discrete mode, the growth rate of which is equal 
to zero at small ar, and becomes negative with increasing 
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Figure 1. The dependence of the growth rate Im (tl>) vs. dimen- 
sionless angular momentum dispersion aj of the mode / = 3 
for models n = 2 (diamonds) and n = 3 (circles). Dashed lines 
show the asymptotic behavior obtained using spoke approxima- 
tion equation J2~25t : hn{p/a T ) = (2/7r 2 )ImA = 0.189 and 0.532 
for n = 2 and 3 respectively (exact solution for aj> = 0.003 gives 
0.185 and 0.529). 



qt- This is just that weakly decaying mode with the fre- 
quency u) 2 w 2 Ci I (I + l)/7r 2 (at oit <C 1) which was men- 
tioned in Sec. 2.2.2. The dependence of the frequency on I 
for this mode was one of the tests for numerical code of the 
integral equation solver. Another test was detecting the zero 
lopsided mode I — 1 mentioned in Sec. 2.2.1. 

The third test was the evaluation of Q (qt) depen- 
dence in the spoke approximation limit. Assuming that 
uj = 2Aqt/vt 2 , and using K ,{a, a') ~ 1, (f> s (a) ~ 1, and 
f(a) w 2a /-k 2 in (|2.2p . one can obtain the equation for the 
/ = 3 mode 



oc 

/ 



dz ( 



A 2 



+ 



15 



A 2 -9z 



= 0(a 2 T ). (2.25) 



By setting the r.h.s. to zero, one obtains an unstable mode 
for each n: A = 2.243 + 0.189i for n = 2, and A = 2.592 + 
0.532i for n = 3. The same values obtained from solution 
of the exact integral equation (12. T[) for ar = 0.003 are A = 
2.240 + 0.185i and A = 2.588 + 0.529i, correspondingly. 



2.3.2 The Heaviside model 

The simplest non-monotonic model that allows to progress 
rather far by analytical methods is the model with a piece- 
wise constant distribution over momentum (|2.23|l . In other 
words, we assume the DF to be non-zero only within the 
interval hi < a < h 2 , where it is taken constant. 

When studying stability of discontinuous distributions 
such as (|2.23jl , it is more convenient to start with the integral 
equation in the form (|2.2[l . Substituting (|2.23[l into Eq. (|2.2[). 
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one obtains 



<Ma) 



4C;efii 



h 2 2 -hj sil u 
xK%(a,hi)<f> 3 ,(hi) 



n pr (/ii) hi 



uj 2 - s'^ni^hi) 

OJ 2 - S ' 2 n2 r (/ l2 ) X 



xK.^ 3 ,(a,h 2 )^ s '(h 2 )\. (2.26) 



Let us turn again to the natural slow scale of frequencies 
according (|2.6[) and then substitute in (|2.26[) particular val- 
ues a — hi and a. = h%. For brevity sake, the following 
designations are used: v\ = v{hi), v% = v{h 2 ). We have 



<]>s(hi) = - 
x 4> B >(hi) 

(j>s(h 2 ) 



iCi 
hi - h\ s , 

V 2 hi 



v\ hi 



K,M,(hi,hi) x 



UJ* — S' " V 



/ 2 , ,2 s s 



LU 2 — s' 2 



(2.27) 



4C; 



h 2 -h 2 

z 1 s' = s 



Vihi m . . 

-2 TT^K- ss>{h 2 ,hi) x 

OJ 2 — s' 2 v{ sa 



(2.28) 



This set of equations relative to <j> s (hi) and (f> s (h 2 ), (s = 
1, 2, [| (Z + 1)]) can be reduced to the standard linear set. 
Introducing new unknown functions 

vi hi 



2 — c2 7/ 2 



U) z — S^V 

one obtains 



3 ( ftl )' Ys = 'Ti 2 ~~ 2 . . 2 "MM, 



f-2 2 2\ v . ^ Vlh 

(LU — S Vi ) X s = — 4Cl 



h E s ' 2d >'-- 



hi - hi s , La 
x [fC (l J sl (hi,hi) X s , - fC { ^ s ,(hi,h 2 )Y a ^ , (2.29) 



i>2 h 2 
hj^h 2 



x [K, {l X,{h 2 ,hi)X s , -K%(h 2 ,h 2 )Y s ,] . (2.30) 

The precession rates in these equations can be expressed 
through the complete elliptical integrals K and E: 



„ 9r hi 1 - q 2 Q(q) 



r> n h 2 Q(q)-q 
ei 1-q 2 



where Ci = 4/(3tt 2 ), q = e 2 /ei, ei = (1 - h\) l/2 , e 2 = 
(1 — /il) 1 ^ 2 , and the function Q(g) is 

Qil) = lh 2 [(l + q 2 )E(q)-(l-q 2 )K(q)]. 
z q 

In the limit h 2 — > hi, the frequencies vi and v 2 are coin- 
cident: vi — v 2 = (2/n 2 ) (h/e), where h = hi = h 2 , and 
e = ei = e2- Note that a determinant of the set of equations 
^T2^i and (^3U|) has a rank 2 [§ (I + 1)]. Particularly, for the 
mode 1 = 1, the rank is equal to 2. Roots of the determi- 
nant are calculated numerically. The difference h 2 — hi has a 
meaning of dispersion, i.e., it is analogous to the parameter 
qt in our models with smooth distributions. 



A simple analytical task is to ascertain that lo 2 = is 
the eigenvalue of the determinant for 1 = 1. We have for 
(«,«') 

/Cg ) (a,a') = e<Q(«;), (2.31) 

where k = e< /e> , e< = min (e, e'), e> = max (e, e'), e = 
(1 - a 2 ) 1 / 2 , e' = (1 - a 12 ) 1 ' 2 . In particular, 



IC^ii' (hi,hi) = ei, A^^fe, fe) = e 2 , 



^(fti,/i 2 ) =/C^ ) (fe,/ii) = e 2 Q(q), 



(2.32) 



q = e 2 /ei. (2.33) 

Setting lD 2 = in the determinant of the set (|2.29|l and 
(|2.30[) . and using the expressions for the elements of the 
kernel (|2.32p and (|2,33[l , we can show that it is equal to zero 
identically. This just means the occurrence of a zero mode 
in the spectrum. Another root uu 2 for I — 1 is positive for 
any values of hi and h 2 , which agrees with our previous 
conclusion (Paper I) that the instability is absent for the 
mode I = 1. 

It is useful to derive equations in (|2.29|l and (|2.30[l in 
the spoke limit, when the distribution is localized in a region 
of small a. This means that we suggest hi <C 1, h 2 -C 1, 
hi < h 2 , set in §FI5\ <f> s (a) = (-l) s , K® s , = (-l) s+s \ and 
finally obtain 



1 = 



4Ci 



h\ 



hi 



J2s' 2 Df 



vi hi 



v 2 h 2 



(2.34) 



For the precession frequencies vi and v 2 , one has in this limit 
vi = (2/n 2 ) hi, i/ 2 = (2/n 2 ) h 2 . Introducing uj = 4\ tv 2 lu, let 
us write down, e.g., the characteristic equation for the mode 



I = 3. In this case D\ 



3_ 



Dl 



r, hence 



i = - 



3tt 2 C 3 



h\ 



hi 



hi 



hi 



+ - 



15/i? 



9h 2 



14 



9h 2 



(2.35) 



12 — hi <§C 1, the role of "self- 
gravity" may be made sufficiently large in spite of small 
parameter C3. This may give the oscillating instability un- 
der certain conditions for hi and /12. The limiting solutions 
serves a test for the model with arbitrary hi and /12. 

The results for the modes I — 3, I = 4 and / = 5 are 
presented in Fig. [2)3-/. They show the boundaries of insta- 
bility domains on the plane (J12 — hi, hi). Left panels show 
the results of computations from the exact set of equations 
()2.29|) and (|2.30[) ; right panels - the results from the spoke 
equation for this model (|2.34|l . It is seen that in the region 
hi <g; 1, /12 — hi <C 1, the results obtained from the spoke 
equation and those from the exact equations do coincide. 
Location of the growth rate maxima at Figs.[2]a, lc, le, as 
well as the values of growth rates are practically the same[f| 
We conclude that for the non-monotonic DF the instability 
is insensitive to a number I of the mode. This is a charac- 
teristic feature for the loss-cone instability. Recall that in 

3 Two additional instability domains for 1 = 5 are explained by 
more complicated structure of the characteristic equation for this 
mode, compared to the modes 1 = 3 and 1 = 4. 
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h 2 -h, 11,-11, 

Figure 2. The stability boundaries and the isolines 10 2 Im (ui) 
on the plane (hi — hi, hi) for the Heaviside model (left: exact 
calculations, right: spoke approximation calculations): a and b - 
for the mode I = 3; c and d - for the mode I = 4, e and / — for 
mode Z = 5. The spoke approximation is reliable in the lower left 
corner of the domain. It is also seen that for I = 3, the growth 
rate sharply decreases when the ratio of the difference hi — h\ 
(the analog of the dispersion ctj< for models with smooth DFs) to 
the size of the loss cone, hi, becomes greater than 2.2. 



models with monotonic distributions, the destabilizing term 
quickly decreases with increasing spherical number I of the 
mode. 



2.3.3 The log-exp model 

In some numerical models (see, e.g., Cohn, Kulsrud, 1978, 
Berczik et. al. 2005, Spurzem et. al. 2005) the initially 
isotropic distribution transforms under action of a massive 
black hole into one monotonically increasing with angular 
momentum, f(a) oc \n(a/h), where h defines the minimum 
angular momentum of a star which is not absorbed by the 
black hole. In this section, we consider the stability of DF 




a T 



Figure 3. The mode 1 = 3 for the log — exp model. Left: isolines 
10 2 ImiI> on the plane (otT,h). Right: the ratio of the growth 
rate Im (ill) to the real part of the frequency Red), for values 
h = 0.1, 0.2, 0.3. 

((11241) . For it, one finds 

l 

i f daa ( ' o?\ l r ., 2/ 2s, 

h 

+ Ei(-l/c4) + ln(ft 2 ) e - 1/ ^] , 

where Ei(z) is the exponential integral. Density profile is 
obtained from (|2.12[) : 

Much as in the power -exp model considered above, the cal- 
culations of the precession rate should be performed numer- 
ically. 

The qualitative pattern of the spectrum for this model is 
similar to that of the power - exp model: when the dispersion 
is not too large, two discrete modes occur, one of which being 
unstable. 

In Fig. |3 (left panel), the isolines 10 2 Imu on the plane 
of parameters (a, h) for the mode I = 3 are presented. Com- 
parison of the figures [2] and [3] shows qualitative coincidence 
of growth rates behavior on the dispersion of DF and the 
size of the loss cone in this and Heaviside models. The right 
panel shows the ratio of imaginary part to real part of Q vs. 
dimensionless angular momentum dispersion olt for several 
values of loss cone size parameter h. 

2.4 Stable models 

Instability of spherical clusters around massive black holes 
was first studied by Tremaine (2005). He considered distri- 
butions of the form 

/(/i,J a )DcJi 6 lnf^-J (2.36) 

in the domain 7 m i n ^ Ii ^ /max, I2 > hl\ (and zero outside 
this domain); I r , I2 = L are the action variables, I\ = I r +L; 
b and h are the real parameters. In the distribution, the loss 
cone is empty for dimensionless angular momentum a < h. 
Tremaine studied the most large-scale perturbations with 
the spherical indices I = 1 and I = 2. 



© 2006 RAS, MNRAS 000,ITltT5l 



Effect of Angular Momentum Distribution on Loss-Cone Instability 9 



In this section we consider two monoenergetic models. 
The first one is 



/(a)=JVln(^), h 



< a < 1, 



(2.37) 



where h < 1 characterize size of the loss cone, N is the 
normalization constant. Dependence of distributions (|2.36|) 
and with f(L) from (|2.37[) on the angular momentum 
is identical. Just this dependence is crucial for stability or 
instability of each specific distribution. Stability of distribu- 
tion 1)2. 37|) for arbitrary values of I is proved in Sec. 2.4.1. 

Another distribution (Sec. 2.4.2) is the simplest mono- 
tonic Heaviside model in a form of the step-like DF, 

/(Q) = T~tf H(a-h) H(l -a). (2.38) 

The factor H(l — a) is added to reflect that the DF domain 
is bounded by circular orbits. 

In Sec. 2.4.3, we prove the stability of spherical systems 
(in the field of a central massive body), all orbits of which 
are circular. 




Figure 4. Spectrum for the log model with h = 0.1, 0.2, 0.3 for 
1 = 1. In all cases the eigenfrequencies are neutral and consist of 
one zero mode (circled point) and continuous part of spectrum 
(points run into a line). To save room the spectra are shown in 
a single plot, but vertically separated from one another (h grows 
from bottom to top). 



2.4-1 The log model 



Distribution in the form (|2.37|l allows to calculate the den- 
sity p(r) explicitly. Using (12.12[) one obtains: 



p(r)=N 



d(a 2 



hi 



(f). 



(2.39) 



where the normalization constant satisfies the relation: 
i 

A 

h 



From the condition h 2 ^ aL(r) = (ir/R) (1-r/R), it 



follows that 

i? min = ii^l-TT - ~h?), 



Rrr 



= + Vi -ft 2 )- 



In presence of the loss cone, h > 0, the radius of the system 
J?max is less than the apocenter radius for a radial orbit, R, 
since stars with low angular momentum, a < h, are absorbed 
by the black hole. Integration (|2.39|l gives for the density 



p(r) = ~ [vV(-R-r) x 



x In 



y/r (R-r) + \/(r - R mitl )(Rm*x - r) 



J -Rmin 



Rn 



- \f(r - -RminX-Rmax - r)J . 

As is seen the density vanishes smoothly at the boundaries 
of the spherical layer r = i? m i n and r = i? ma x. The expres- 
sions for the precession velocity are defined by the formulas 
(pT0l) - (pTT1) . 

In Fig. [4] the frequency spectrum, w, of the spherical 
harmonic I = 1 for the log -model is presented for differ- 
ent values of the parameter h. All calculations detect zero 
modes. 

For other values of I (we considered I = 2, 3), the dis- 
crete modes are absent in the frequency spectra for all h. 



The spectra are continuous and lie at the region of real and 
positive values of uj 2 . We conclude that the log models are 
turned out to be stable. 



2.4-2 The monotonic Heaviside model 

Following the procedure described in Sec. 2.3.2 for the unsta- 
ble Heaviside model, one can derive the following equation 
for the distribution function (|2.38[) : 



h 2 s'=s 



K { ? sl {h,h)<i> s ,{h), 
(2.40) 



where A = uj/v{h), v(K) = 8h/[3 iv 2 e(h)]. It is easy to 
see that for I = 1 there is a zero mode only, since 
IC\i(h,h) = e(h) = (1 — h' 2 ) 1 / 2 . Introducing new variables 

X s = scj) s (h)jDf /(A 2 - s 2 ) one can reduce (|2.40j) to a 
standard linear set 



X 2 X s = s 2 X s - WAQX'" 



(2.41) 



where the matrix 



Ci KUh,h) 

Ci e(h) 



D? Df 



is Hermitian. So, the eigenfrequencies A 2 are real. Here again 
we have a competition of opposite factors, expressed by the 
first and the second terms in the r.h.s of (|2.4ip . To conclude 
whether the instability occurs, we must find numerically ze- 
ros of the determinant ||(s 2 — A 2 )5 aa / — -£^/|| = 0, as a 
function of A 2 for I ^ 2. A rank of the determinant is equal 

to [ia+i)j. 

The results are represented in Figs. [5]and[6] This figures 
show the dependence of \ 2 (h), j = 1, 2, [| (I + 1)] as a 
function of the loss cone size h. The instability is clearly 
absent. Each eigenvalue A 2 (/i) has only a weak dependence 
on h and approximately equals to j 2 . The least stable mode 
is I = 3 (see Fig. |6j, but it is still far from instability. 
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Figure 5. The dependence of the eigenfrequencies squared, 
A|(/i), for I = 1,2,3,4. 
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Figure 6. The dependence of the smallest (for given I) eigenfre- 
quencies squared, A 2 (/i), for 1 = 1, 3, 5, and 7. 



2.4-3 Model with circular orbits 

Let us consider the simplest monotonic model, in which all 
orbits are circular. In this section we do not assume a distri- 
bution to be monoenergetic, otherwise the density distribu- 
tion would be degenerated to a thin spherical layer. Let us 
assume that the DF is F(E, L) — A S [ L — L c i TC (E)] , where 
A = const is normalization factor and E changes in some 
range AE. In terms of radial and transverse velocities v r 
uE, the DF is: 



and v± = ^/i 



F(v r ,v ± ,r) = -^^6(v r )5[v ± 
2iv vo(r) 



vo{r)], Vo(r) = rfi(r), 



where O(r) is the angular velocity of a star on the cir- 
cular orbits. This velocity is determined by a balance of 
the centrifugal force and a sum of the gravitational forces 
from the central body and from the spherical cluster: Q 2 = 
&o(r) + (l/r)d<S> G (r)/dr, fi§(r) = GM c /r 3 . Here we also 
suggest that Qq » r -1 d<&G(r)/dr. 

In this approximation, the orbits are near-Keplerian, 



wo = 47rGp (r) = $ G + - <E> G , 
r 

1 



fi = Q + 



2rQ, 



(2.42) 



2i 2o V r 



For the precession rate we have (see also Tremaine, 2001) 
fi pr = Q, - k = -(l/2fio) [*g + (2/r) <E> G ] , or, taking into 
account ((2^421 , fl pl = -|u§/n . Since e = M G /M C < 1, 
one has the following scalings: S7 pr ~ ef2o, w§ ~ eO 2 ,, and 
for slow modes w ~ Q pr ~ efio- 

We start from the equation derived by Pal'chik et. al. 
(1970), for the models with circular orbits (the equation 
can also be found in the monograph by Fridman and Poly- 
achenko (1984)Q It has a form: 

J. r 2 A, (r, w)^ - B, (r, w ) » (r) = 0, (2.43) 

where xK 7 ") i s a radial part of the potential perturbation 
$(r) oc Xi( r ) ^; m (^! ex P( — iwi), coefficients Ai(r,ui) and 
Bi(r,u>) are 



A 8 



[oj— (sfl — k)][uj — (sfi+Ac)] 



(2.44) 



B,(r,u;) = Z(Z + l) + E °? 



i 2 
+ W 



9 
Wo 



2sf7 



r (ui — - sfl) [w — (sf2 — t)][w — (sfi + Ac)] 
4sf2(a; - sfi) + s 2 [(w - sfi) 2 + 4fi 2 - k 2 
(w - sfi) 2 [w - (sfi - k)][lo - (sQ + At)] 

(i + s + l)0-s) 



-+ 



(w - sfi) (w - sfi - 2fi) 



(2.45) 



Now we need to distinguish between even and odd val- 
ues of Z, since I and s should be of the same parity (i.e. both 
even or both odd). 

For even I, the dominating contributions is expected 
from s — and s = — 2. However, one can see that for even 
I, the contributions from s = and s = — 2 cancel each 
other. Indeed, setting cjq = — 2ilo^pr one has: 

A; = 1, B, = J(f+1) + D? J(J + 1) - A 2 (/-I) 
After taking into account the relation 

Df = A°{J(J + l)/[(J-l)(J + 2)]}, 



one obtains: ^4; = 1, B; = Z(Z + 1). So the equation (|2.43|l 
is reduced to the trivial relation A\i = 0, which means the 
absence of the slow density perturbations. 

For odd I the terms s = ±1 give the main contribution to 



4 Note that both in the monograph and in the original paper, 
the form of equation does not allow to include the external grav- 
itational field from a halo or a central body. We have slightly 
changed the equation to make it possible. 



© 2006 RAS, MNRAS OOO.riHTsl 



Effect of Angular Momentum Distribution on Loss- Cone Instability 11 



the sum, while other terms (|s| 7^ 1) are beyond the accuracy 
of the slow mode equation. Thus, one has: 



Ai = 1 + 2 Dj 



fin 



On 



dr \r - J7p r 



(2.46) 



To study this case we transform the differential equation 
(|2.43|) with Ai and Bi from (|2.46[1 to an integral equation. 
Eq. (|2.43[) can be represented in the form of the Poisson 
equation 



A X i(r) =4ivGpi{r), 
with the perturbed density 



(2.47) 



(2.48) 



where S{tu 2 ,r) = 2 n£ r /(w 2 - Q% T ). The solution of Eq. 
(12.471) in the integral form is: 



Xi(r) 



4ttG 
21 + 1 



r' 2 dr' pi{r')Ti(r,r') 



(2.49) 



where the kernel 



*i(r, O = ^lH{r- r>) + ^ H(r> - r), (2.50) 
or, substituting (|2.48[) into (12 .49 p . and integrating by parts, 

Xi(r) = - A / dr'S{u 2 ,r') A [/ 2 Xi (Ol x 



dfi(r,r') 2 , 
+ 7^) 



(2.51) 



Applying the operator "P(r) = d/dr + 2/r to both parts of 
Eq. (p3Tj) and denoting *;(r) = V(r)xi(r) = dxi(r)/dr + 
(2/r) Xi( r ) we obtain an integral equatioqfl 



*i(r) = - 



_2L 

2/ + 1 



r' 2 dr'SV 2 , r') ft; (r, r') * ; (r), (2.52) 



with the new symmetrical kernel IZi (r, r') = 

(d/dr + 2/r) (d/dr' + 2/r') Ti(r,r'). Introducing the 

new function Zi = r Q pr /(uj 2 — fi pr ) ^i, we obtain the 
required integral equation 



K - fi pr (r) 2 ] Z,(r) 



2D} 
21 + 1 



dr' fi pr (r) Q pr (r') fci(r, r') Z,(r'), (2.53) 



with the kernel ICi (r, r') = rr' 1Zi(r, r'). 

Since the kernel defines a self-adjoint integral operator, 
all eigenfrequencies uj 2 should be real. To determine whether 



5 The integral equation H2.52H can also be derived from the 
general "slow" integral equation, by considering of circular orbit 
limit. However, that derivation is much more cumbersome than 
one given here. 



negative values of uj 2 are possible, let us write out the kernel 
K,i(r,r') explicitly: 

Ki(r,r') = -{I + 2) (I- l)Fi(r,r') + (21 + 1) 5(r - r'). 

(2.54) 

It contains two contributions: the first is negative, the second 
is positive. Substituting (|2,54[) into (|2,53|) . one finds 

u 2 Z l {r) = fi. 2 r (r)(l - 2D}) Z,(r) +2 D\ ^-M^ x 
x / dr' n pr (r)n pI (r') fi(r,r) Zi(r'). (2.55) 



For / = 1 one can see that uj 2 — satisfies this equation 
since D 1 = 5. However, the most interesting fact consists in 
stability of all higher modes, I 3. Indeed, since l-2D\ > 
for I 3, and the integral operator in the r.h.s. is positively 
defined, we conclude that all eigenvalues, to 2 , are positive. 
Consequently, the instability is absent in the limit of circular 
orbits. The result is universal and does not depend on a 
particular choice of the model Q pr . 



3 THIN DISK SYSTEMS 

3.1 Slow mode Integral equation for 
monoenergetic disk models 

In this section we shall consider the monoenergetic distribu- 
tions of (|1.1[) type assuming the function f(a) to be even, 
f(a) — /(—a). The function F(E,L) is normalized as fol- 
lows: 



J F dT = (2tt) 2 y" 



^circ(S) 



dLF{E, L), 



-L cizc (E) 



which gives the normalization constant 

A= Mo 
ir 2 R 2 



(3.1) 



(3.2) 



provided that J*_ a f(a) da — 1. 

The integral equation for slow modes (Paper I) can be 
represented in the form: 



C m fda'df/da' 



TV L ' J LO — V 

-1 



-—— Km (a, a) 4>{a), (3.3) 
(a') 



or, using the evenness of <f)( a ) which stems from evenness 
of f(a), oddness of fi pr (a), and symmetry properties of the 
kernel, K, m (a, a') = K, m (— a, a'), 



0(a) 



2 CV? 



7T" J UJ 




v{ct) df 



____£ m(Q ,a')0(a')da' 



(3.4) 



where uj and v{a) are the dimensionless pattern speed and 
the dimensionless precession rate: 

fip / *\ fi P r(c\) , . 

w=— — , i/(a) = — — . (3.5) 
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Changing the unknown function, integral equation (|3.4p 
takes the form of linear eigenvalue problem: 

1 

[lu 2 ~ v 2 {a)] ip(a) = 2 ^ m I v{a) — K m {a,a) ip(a) da'. 



(3.6) 



The kernel functions for thin disks K, m {a>, a') can be trans- 
formed from the corresponding expression in Paper I to a 
suitable form as follows: 

7T 7T 

fC m (a,a) = —J— / drr cosm( / dr' r cosm(' jT m (r,r'), 
Gm J J 



(3.7) 

where dependence of r and anomaly £ on r and e are the 
same as in the spherical case (|2.8|l , but the function T m (x, y) 
is 



Tremaine 2005 and Paper I) 

n v 8 1 



2 d$ G , 
r — ; — cos (, dQ. 



dr 



(3.14) 



but the relation between the potential and the surface den- 
sity is much more complicated^] (see, e.g., Tremaine, 2001) 



$ G (r) 



AG 

-1/2 



(r') 1/2 <T (r')dr' [q 1/2 K(q)]. (3.15) 



Using (|3.12p - (|3.15[) one obtains a suitable expression for 
the scaled precession rate v(a) (|3.5|) in the integral form 

l 

v{a) = a j dot f(a) Q{a, a), (3.16) 
o 

where Q(a,a') is a universal function (i.e. does not depend 
on form of the distribution): 



cos mO d9 



n yj ' x 2 + y 2 — 2xy cos 6 



(3.8) 



As before, the kernel /C m (a,a') is normalized to unity: 
IC m (0,0) = 1, which means C m equal to 



dx 



1-x 



dy 



T m {x,y). (3.9) 



The latter formula immediately follows from (|3.7[) if one 
reminds that for radial orbits £ = n, cosm( = (—1)'", 
dr = dx [x (1 — x)]" 1 / 2 . For the lowest azimuthal numbers, 
functions J- m (x,y) can be expressed through elliptical inte- 
grals of the first and the second kind K(q) and E(g): 



Pi{x,y) 



4 K(q)-E{q) 



r> 



(3.10) 



^2{x,y) = - — 
3r> 



^ + l)K( g )- 2 (i+l)E( 7 ) 



(3.11) 



where r> = max(x,y), r< = min(a;,y), q — r < /r > . Using 
(|3.10l) and (|3.11[l one can obtain numerically C\ = 10.88, 
C 2 = 7.45. 



For the surface density we have: 



ao(r) = I JdE J 



F(E, L) dL 



2M G 



Q(a,o') 



r' dr' 



x dr ■ 



% 5 r / n ^/(r'-r^J^-r') 
E(k) K(k) 



2r - a 2 



■ y (r— r m in)(r max — r) 



r r' + r 



(3.17) 



Here k = 2 V r r'/ (r + r'). The integral of the first term is 
understood in the principle value sense. Using the same trick 
as in Sec. 2, one can change to new integrating variables r 
and r', where r = | (1 — e cosr) and r' = i (1 — e' cost'). 
Then, for Q(a,a') one obtains 



Q( 



a, a = 



2tt 3 ( 



e — cos r 
1 — e cos r 



dr' (1 — e' cos r ) x 



E(k) K(k) 



r' + r 



(3.18) 



where 



-iWM \/2S + 



a max (r) 



2GM C 



So(r), 



So(r)= | 



/(a) da 



-a max (r) 



(r) — a 2 



TT 2 R 2 

(3.12) 
(3.13) 



3.2 Variational principle and sufficient condition 
for instability of m — 1 mode 

As we see from Sec. 2, for spherical systems with the mono- 
tonic DF, the variational principle takes place. Besides, for 
/ = 1 and the empty loss cone, a zero frequency solution 
exists which stands for a sphere displacement from the mas- 
sive center; all other eigenmodes being stable. A thin disk 
is completely different. The displacement is no longer an 
eigenmode. Moreover, models with analogous distributions 
are turn out to be unstable. Let us prove the instability of 
lopsided m = 1 mode provided that 



(a) = u{a)df/da < 0. 



(3.19) 



and «L(r) = ^maxM / L 2 ilc = 4(r/R)(l - r/R) as for 
spheres. 

The relation between the precession rate and the po- 
tential &c(r) is the same as in the spherical systems, (see 



Note that spherical models with the analogous condition are 
stable. 

Disks with even DFs satisfying condition (|3.19[l obey 
the variational principle, which means the eigenfrequencies 
squared are real. So one can formulate a sufficient condition 

6 For that reason the precession in the near-Keplerian disk is not 
always retrograde. 
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of instability for m — 1 azimuthal perturbations as follows: 
If the loss cone is empty (/(0) = 0), the DF is monotonically 
increasing df/d\a\ > 0, and the precession is retrograde for 
all values of angular momentum (v(a)/a < 0), then m = 1 
perturbations are unstable. 

To prove the statement we shall use integral equation 
(13.41) . in which ui = iT with T > is assumed: 

M(T)<j>(a) = 0, (3.20) 

where operator -M(F) is 

i 

M(T)J>(a) = 4,(a) + ^r / r2 f Q 2 ,M a,a')<j>{a')da'. 
T J 1 + f (a J 

(3.21) 

Now we consider another eigenvalue problem 

M{T)<j)(a) = X(T)cj>(a). (3.22) 

Eigenvalues F of the problem (|3.20[) correspond to eigenval- 
ues A (F) = of the problem (|3.22p . Let us define an inner 
product as {X, Y) = J* daX* (a)W{T,a)Y{a), where the 
weight function W(T,a) = -g(a)/[T 2 + v 2 (a)} > 0. Opera- 
tor M has the following properties. 

1. A4 is Hermitian, i.e. {ip(a),A4(f)(a)) = 
(Mi>(a),<f>(a)) = {<t>(a),MiP(a))*. 

2. M is continuous, when Y 0. One might think that 
the first term (j>(a) in the r.h.s. of (|3.21[) breaks down the con- 
tinuity, which in turn means that system of proper functions 
is incomplete. However, this is not the case, since 4>{a) can 
be absorbed by introducing new eigenvalue A(F) = A(r) — 1 
in (I3.22p . Since /(«) is even, for smooth DF we have /(0) = 
/'(0) = 0, /"(0) > 0. This condition guarantee the weight 
function W(a) to be finite even for T = 0, despite v — 0(a) 
at a — > 0. 

3. A4 is positive definite at sufficiently large F. This 
is evident, since W(a) > 0, and the second term in (|3.2ip 
becomes small at large F. 

;,From the first two properties it follows that for fixed 
T ^ eigenvalues A„(r), n = 1, 2, 3, ... of A4(T) are real, 
and the system of proper functions is complete. The third 
property means that at large F all eigenvalues A n (F) are 
positive. 

If we find a test function 4>t(To, a), for which the scalar 
product (0t(Fo, a), A4(To) 4>t{To, a)) is negative, it mean 
A4(To) is not positive definite for the given Fo. So, at least 
one eigenvalue must be negative: A m i n (ro) < 0. This mini- 
mal eigenvalue A m in(r) increases with V and becomes posi- 
tive, as all other A n (F). We conclude that there must be a 
value of r, To < r < oo for which A m i n (r) = 0. This value 
is an eigenvalue for f|3.21 f) , which means the existence of the 
eigenmode describing aperiodic instability with growth rate 

r. 

For the test function 0t(Fo, a) one can take a displace- 
ment of the disk from the center, which is similar to the 
sphere displacement (|2.15[) and correspond to the lopsided 
perturbation m — 1: t (Fo,a) = (e/a) v(a), and Fo = 0. 
One can show that 

{4> t {T ,a),M{T )ct>t(T ,a)) < 0. (3.23) 



Let the l.h.s. be —P, or, explicitly 

n f , df(a) re(a)l 2 . . 2Ci } J e(a) 

P= I da , ' — — via) H / da— — x 

J da |_ a J tt J a 

o o 

df(a) f,,e(a') df(a') 
o 

After some lengthy manipulations using (|3,7[) , (|3.10|) , (|3. 18[) , 
and condition /(0) = 0, one can show that P is positive, so 
inequality (|3.23|) is fulfilled. 

Tremaine (2005) has also obtained a sufficient condition 
for a lopsided mode in the symmetrical disk using Good- 
man's (1988) criterion. His condition, however, differs from 
ours. Namely: If the loss cone is empty, F(E, L = 0) = 0, 
and d^a{r) / dr > throughout the radial range containing 
most of the disk mass, then disk is unstable with respect 
to m = 1 perturbations. This formulation does not use the 
requirement for precession to be retrograde and the DF to 
be monotonically increasing, although monotonic increase of 
F(E,0) = [8F(E,L)/dL\ L=0 = 0, [d 2 F(E, L)/8L 2 ] L=0 > 
is implied, at least for small angular momentum. Thus, the 
comparison between spherical and disk case can hardly be 
made, unless the conditions of stability is formulated in sim- 
ilar terms. To perform the comparison, we give our own cri- 
terion that follows directly from the integral equation. 

It needs to be emphasized that the sufficient condition 
by Tremaine (2005) is different. Lack of the condition for the 
sign of precession possibly means that his criterion include 
two type of instability simultaneously - the radial orbit in- 
stability arising in disks with prograde precession, and the 
loss cone instability which require retrograde precession. 

For disks composed of near-radial orbits, Tremaine's 
condition gives the result obtained in Paper I: a disk with 
symmetrical distribution f(a) obeying the conditions /(0) = 
/'(0) = 0, /"(0) > is unstable if the precession is retro- 
grade. In turn, precession of near-radial orbits is retrograde 
if d$ G {r)/dr > 0. 

3.3 Numerical results 

To support the mathematical rationale given above and to 
provide a basis for possible simulations, it is useful to obtain 
eigenfrequencies of unstable modes for particular models. 
Here we consider the power -exp model with symmetrical 
distribution 

/(a) =7v4-exp(-a 2 A4), (3.24) 
where the normalization constant is 

TV -1 = 2 J x 2 exp(— x 2 ) dx. 
o 

For qt <C 1 the constant TV = Distributions become 

monotonic in the interval [0, 1] when ar ^ 1. Note, that 
when cut 3> 1 DF is simply /(a) = § a 2 in the interval 
[— 1, 1] and doesn't depend on aj. 

Evaluation of integral equation (|3.6|) requires prelim- 
inary calculations of the kernel function JC m (a,a') using 
(|3.7[) . and the scaled precession rate v(a) using (|3.16[) and 
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Figure 7. The dependence of 7 (growth rate divided by az- 
imuthal number m) vs. dimensionless angular momentum disper- 
sion cut of the initial DF for azimuthal numbers m = 1 (crosses), 
m = 2 (circles), and m = 3 (squares). 

(|3.18[) . For brevity, we skip the details here, just noting that 
calculation of function Q(a,a') is turn out to be rather dif- 
ficult numerical task. The calculations show that for the 
model (I3.24p the precession rate v(a) is retrograde for all a, 
i.e. v{a)/a < 0. 

The results for values of azimuthal number m = 1, 2, 3 
are collected in Fig. [7] Since the initial distribution is sym- 
metric, real parts of the eigenvalues u> are equal to zero. 
Hence in Fig.[7]we show the imaginary parts 7 = Imw, which 
are the growth rates of the unstable modes divided by az- 
imuthal number m, vs. dimensionless angular momentum 
dispersion «t. One can see that instability exists for all ar 
and never becomes saturated. Moreover, it is easy to ob- 
tain the asymptotic values 7 for different m at cyt — > 00: 
0.289, 0.108, and 0.026 for m = 1,2,3 correspondingly. 
For small angular momentum, growth rates increase linearly 
with QfT, such as f/ar are equal to 0.454, 0.463, and 0.481 
for m = 1, 2, 3 correspondingly. 



4 DISCUSSION 

We have studied the stability of the spherically-symmetric 
and thin disk stellar clusters around a massive black hole. 
We conclude that stability properties of spherical clusters 
depend crucially on monotonity of initial distribution func- 
tions, while thin disk clusters are almost always unstable. 

If the initial distribution of the spherical cluster is 
monotonic, the cluster is most likely to be stable. This con- 
clusion was first made in Tremaine (2005), where stability of 
I — 1 mode was generally proved, and I = 2 was tested nu- 
merically. We confirm this conclusion by considering a num- 
ber of monotonic distributions for modes with arbitrary 
Besides, we have checked distributions obtained from mono- 
tonic ones by making them vanish quickly but smoothly at 
circular orbits. These models were also stable. However, a 
general proof of stability for any monotonic distributions 
was not yet found. 

Spherical clusters with the non-monotonic DFs should 



be generally affected by the gravitational loss-cone instabil- 
ity. The instability was first found in our Paper I using a 
simplification of systems with near-radial orbits. In the Sec. 
2 we show that this instability is due to just non-monotony 
of distributions over angular momentum, the orbits may not 
necessary be near-radial. 

In our opinion, both monotonic and non-monotonic dis- 
tributions are important for possible applications to real 
stellar clusters around black holes. The DFs monotonically 
increasing from the loss cone radius up to circular orbits 
are formed naturally due to two-body collisions of stars. 
It follows from numerical experiments (see, e.g., Cohn and 
Kulsrud, 1978), which predict establishment of such distri- 
butions after a characteristic time for collisional relaxation. 
These distributions may be approximated by the formula 
F(x\n(L/L 

min ) ■ 

Such a slowly increasing function is, in fact, predeter- 
mined by the boundary conditions imposed in the cited nu- 
merical study and some other investigations. Indeed, the 
vanishing condition at L = L m t n , and the matching con- 
dition to isotropic (Maxwellian) distribution, F = F(E), at 
the boundary E — Bbouad = of the phase space (E, L) 
(boundary separates stars which is gravitationally coupled 
to the black hole from the others) is required. The last 
condition means the asymptotic (when E — > Abound) in- 
dependence of the function F(E, L) on the momentum L. 
So monotonic, or logarithmic, dependence of type of (|2.37[1 
is quite reasonable. 

The non-monotonic distributions are also real. If the 
cluster, is formed, for example, as a result of the collisionless 
collapse (several free fall times), then it remains collisionless 
for a long timescale of collisional relaxation (see, e.g., Mer- 
ritt & Wang, 2005). In principle, the system can have almost 
arbitrary DF both in the energy and in the angular momen- 
tum. During the collapse, a typical non-monotonic distribu- 
tion of stars over the angular momentum, with empty loss 
cone and maximum at some value L = L*, is formed. 

In Paper I we argued that stability properties of such a 
distribution is effectively analogous to one of typical plasma 
distributions of the "beam-like" type. But they can read- 
ily become unstable, as it is well-known in plasma physics 
(and also confirmed by direct stability study of correspond- 
ing stellar systems in Paper I). It is possible (as it is often so 
in plasma) that for the time of collisionless behavior, DF can 
undergo a dramatic change from its initial form. In partic- 
ular, the collective flux of stars into the loss cone caused by 
the instability could, in principle, lead to the formation of a 
considerable part of the black hole. Checking of such possi- 
bilities is the most urgent task for future studies of unstable 
non-monotonic models. 

Since spherically-symmetric models with the monotonic 
DF are apparently stable, but analogous disk systems are 
unstable (see Tremaine 2005 and Sec. 3), a critical flatness of 
ellipsoid models at which the instability begins is expected. 
Study of such systems, as well as systems with more complex 
triaxial ellipsoids can be performed using numerical simula- 
tions. 
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